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Abstract 

We analyze the possible phase diagrams of a simple model for an associating liquid proposed 
previously. Our two-dimensional lattice model combines oreintational ice-like interactions and "Van 
der Waals" interactions which may be repulsive, and in this case represent a penalty for distortion 
of hydrogen bonds in the presence of extra molecules. These interactions can be interpreted in 
terms of two competing distances, but not necessarily soft-core. We present mean-field calculations 
and an exhaustive simulation study for different parameters which represent relative strength of 
the bonding interaction to the energy penalty for its distortion. As this ratio decreases, a smooth 
disappearance of the double criticality occurs. Possible connections to liquid-liquid transitions of 
molecular liquids are suggested. 
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I. INTRODUCTION 



The recently acknowledged possibility of the existence of single component systems which 
display coexistence between two different liquid phases has opened many interesting 

questions regarding the link between the interaction potential and the presence of double 
criticality . 

Network-forming fluids are primary candidates for liquid-liquid transitions. They ex- 
hibit directional, intermolecular attractions that lead to the formation of bonds between 
molecules. As a result, locally structured regions have lower density than unbonded regions. 
These structures at the fluid phase are transient and local but become permanent at lower 
temperatures. 

The case of water is probably the network-forming liquid most intensively studied, due 
to its ubiquity in nature. An unstable liquid-liquid transition ending in a critical point 
was initially proposed to explain the anomalous behavior of network-forming liquids such as 



water jl|| |2lj-p- Closely associated with the possibility of a liquid-liquid phase transition in 
supercooled water is the phenomenon of polyamorphism occurring below the glass transition 
temperature, Tg. Polyamorphism refers to the occurrence of distinct amorphous solid forms. 
There is by now a general consensus that water displays a transition between two different 
amorphous states in the supercooled region of its phase diagram [7]. Experiments carried 
out in water at T 130 K show an abrupt change of volume as a function of pressure which 

n 

indicates the existence of a first-order transition 8]. The two amorphous phases of water 
might be related to two different liquid phases at higher temperatures. At coexistence, 
these two liquid phases determine a first-order transition line ending at a critical point. 
Simulations and experiments predict that the liquid-liquid transition is in an experimentally 
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inaccessible region of the phase-diagram 

Notwithstanding of its confirmation for metastable water, the coexistence of two liquid 
phases was uncovered as a possibility for a few both network-forming and non-bonding liq- 
lations for realistic models for carbon Q], phosphorus [l^ . 



uids. Computer simn 
Q, silica {Si02) 14 1 



germanium 

Q QQ and ...con QQ l^e e'xistenceTf a fl.t-o.de. 

transition between two liquid phases. Recent experiments for phosphorus [s^j^ll or phos- 



phate compounds 



22j |. and for carbon 23] confirm these predictions. 



Substances that are structurally similar to water, such as Si, Ge, Ge02 and silica, can 
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exhibit not only the hquid-hquid transition but also the other thermodynamic and dynamic 
anomalies present in water. Particular attention has been given to silica because its tech- 
nological applications. Many intriguing properties of liquid silica and water occur at low 
temperatures. Examples include negative thermal expansion coefficients and polymorphism 



for both liquids. Compression experiments on silica show a non-trivia 



scopic structure 
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change in the micro- 
28[ , was interpreted 



25| . This behavior, reproduced by simulations 
as an indication of a ffist-order transition. Inspired by these results, new simulations for 
different models for silica were performed giving support to the hypothesis of 

a first-order transition in silica. However, experiments must sill confirm a discontinuous 
.....tionbetwee. these a,„o.phous phase JQ. 

In resume, a full understanding of the effects of the number, spatial orientation, and 
strength of bonds on the global fluid-phase behavior of network fluids is still lacking. 

In order to gain some understanding, a number of simple models have been proposed. 

r~i 

Since the work of Bernal ISOl, the water anomalies have been described in terms of the 
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the presence of an extensive hydrogen bond network which persists in the fluid phase |.3L] . 
In the case of lattice models, the main strategy has been to associate the hydrogen bond 
disorder with bond 



ll|j3^ or site 3^ 3^ Potts states. In the former case coexistence 
between two liquid phases may follow from the presence of an order-disorder transition and 
a density anomaly is introduced ad hoc by the addition to the free energy of a volume 
term proportional to a Potts order parameter. In the second case, it may arise from the 
competition between occupational and Potts variables introduced through a dependency of 
bond strength on local density states. 

A different approach is to represent^drogen bonds through ice vat.abiesHQHQ, 
SO successful in the description of ice j39| entropy, for dense systems. In this case, an order- 
disorder transition is absent. Recently a description based also on ice variables but which 
allows for a low density ordered structure j3| was proposed. Competition between the ffiling 
up of the lattice and the formation of an open four-bonded orientational structure is naturally 
introduced in terms of the ice bonding variables and no ad hoc introduction of density or 
bond strength variations is needed. Our approach bares some resemblance to that of some 
continuous models 



4l| 42]j43|, which, however, lack entropy related to hydrogen distribution 
on bonds. Also, the reduction of phase-space imposed by the lattice allows construction of 
the full phase diagram from simulations, not always possible for continuous models j4l| . In 
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FIG. 1: High density liquid, HDL with density 1 (top) and low density liquid, LDL with density 
3/4 (bottom)on the triangular lattice. The solid lines indicate the hydrogen bonds where the 
arrows differentiate bond donors from bond acceptors. 

our previous publication, we have shown that this model is able to exhibit for a convenient 
set of parameters both density anomalies and the two liquid phases. Here we explore the 
phase space for various values of the orientational energy term. We show that by varying 
the relative strength of the orientational part, we can go at a fixed temperature from two 
coexisting liquid phases as observed in amorphous water to a smooth transition between two 
amorphous structures as might be the case of silica. 

The remainder of the paper goes as follows. In sec. |n]the model is revised, for clarity. In 
sec. mil mean-field results are given. Sec. IIVI has the simulations results. Conclusions end 
this session. 

II. THE MODEL 

Consider a lattice gas on a triangular lattice with sites which may be full or empty. Besides 
the occupational variables, cxj, associated to each particle there are six other variables, 
tI\ pointing to neighboring sites j: four are the usual ice bonding arms, two donor, with 
r*-' = 1, and two acceptor, with r*-' = — 1, while two additional opposite arms are taken as 
inert (non-bonding), r/"' = 0, as illustrated in Fig. 1. Therefore each occupied site is allowed 
to be in one of eighteen possible states. 

Two kinds of interactions are considered: isotropic "van der Waals" and orientational 
hydrogen bonding. An energy —v is attributed to each pair of occupied neighboring sites 
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FIG. 2: Effective potential vs inter-particle distance for u/v=0.5 (dashed line), 1 (full-line) and 2 
(dot-dashed line) . Penalty on non-bonding neighbors is sufficient to stabilize the low density phase 
for u/v=l and 2. The shoulder potential is attractive in the first and repulsive in the second case. 
For smaller penalty u/v=0.5 only the high density phase is stable at zero temperature and the step 
potential is accordingly " smoother" . 

that form a hydrogen bond, while non-bonding pairs have an energy, —v + 2u (for u > 0), 
which makes —2u the energy of a hydrogen bond. The overall model energy is given by 



where o"j = 0, 1 are occupation variables and r*-' = 0, ±1 represent the arm states described 
above. Note that each particle may have six neighbors, but the number of bonds per molecule 
is limited to four. For u/v > 1/2, the "van der Waals" forces become repulsive. As a result, 
each molecule attracts four neighbors, if properly oriented, and repeals the other two. An 
interpretation for this "repulsion" would be that the presence of the two extra neighbors 
distorts the electronic orbitals, thus weakening the hydrogen bonds. 

Inspection of the model properties allows the prediction of two ordered states, as shown 
in Fig. 1. For low chemical potential, the soft core repulsion becomes dominant, p = 0.75, 
and energy "volume" density is given by e = E/V — —3v/2, where V is the number of 
lattice sites. If the chemical potential is high, p — 1, and energy density e = —3v + 2u. At 
zero temperature, the low density liquid (LDL) coexists with the high density liquid (HDL) 



(1) 



5 



at chemical potential /i/f = — 6 + Sw/f , obtained by equating the grand potential density 
(or pressure) associated with each one of these phases. Similarly the coexistence pressure at 
zero temperature is given hy p/v = — 3 + Qu/v. Besides these two liquid states, a gas phase 
is also found and it coexists with the low density liquid at chemical potential [ijv = —2 
and pressure p = 0. The condition for the presence of the two liquid phases is therefore 
ujv > 0.5, i.e., that the "Van der Waals" interactions are repulsive. 

The two ordered structures could be qualitatively associated to different ice phases. Under 
pressure, hydrogen bonds reorganize in more dense phases 

Our model may be interpreted in terms of some sort of average soft-core potential for 
large hydrogen bond energies. The low density phase implies average interparticle distance 
diD = P~LD^ = 2/-\/3, whereas for the high density phase we have duo = Pud — 1- The 
corresponding energies per pair of particles is —v and — f + 2m/3. The hard core is offered by 
the lattice. For ujv > 3/2, the shoulder becomes repulsive, making the potential soft-core. 
Figure 2 illustrates, in a schematic way, the average pair energy Cp/ v dependence on distance 
between particles, for u/v = 0.5, 1, 2. Notice that for u/v = 2 the potential at high density 
phases becomes repulsive. For u/v = 1, 0.5, the shoulder is attractive. 



III. THE MEAN-FIELD ANALYSIS 



Designing a correct mean-field version of the model is not an obvious task, due to the 
special orientational ice-like interactions. We have therefore looked at model properties on 
the Bethe lattice, for which exact relations may be given. Each site of the usual Cayley 
tree (of coordination six) is replaced by a hexagon and {rj, r) variables are attributed to its 
vertices (see Fig. 3), with ri and r defined as in the previous section. This representation 

nn 

is inspired on the Bethe solution for a generalization of the square water model I4fi| 
presented by Izmailian et al. For an occupied site-hexagon, we will have two vertices 

with {ri,T) = (1,1), two with (r/, r) = (1,-1) and another two with (//, r) = (1,0). For an 
empty site-hexagon the six vertices have (r], r) = (0,0). 

The partition function in the grand canonical ensemble is given by: 

S(T, /i) = ISe^^gUh +)9Uh -)9Uh 0) + ^^(0, 0), (2) 
where n is the chemical potential and gfyf{l,+) is the partial partition function (generation 



6 



avo 
^6 




FIG. 3: The Bernal model on the Bethe lattice. Each site of an usual Bethe lattice is placed by a 
hexagon with (ry, r) variables placed on the vertices. Sites of the N generation are linked by dashed 
lines. 

A^) for a central site with (77, /i) = (1, 1). Similarly for (?Ar(l, — ), giy{l,0) and gN{0,0). 

On each lattice line we have two (77, r) pairs of variables connecting two consecutive 
generations of the lattice. The Boltzmann weights for each lattice line are u; = 1, if at least 
one of the connecting sites is empty, and, if both sites are occupied, either uj = e^'", if arm 
variables are of opposite sign, and to = 6^*^^"^"^ otherwise. 

The density of water molecules at the central site of the tree is given by 

18z + x^yj^rj^ 

where z = e^'^ is the activity and 

_ ^?iv(o,o) _ ^7v(o,o) _ gN{o,o) 

9n[1^+) 9n[1^-) 5'iv(l,0) 

Eqn|21 gives us a relation between the water density p and the activity z: 6z = 
lixNyNrN)"^, with 7 = p/(3(l — p)). Using this result the recursion relations at the fixed 
point may be written as 

7(x + y + r) + 1 



X 



^g/3(«-2n) (-J, _^ ^2f3uy _^ ^) _^ ^ ' 

7(x + y + r) + 1 
_ 'y{x + y + r) + l 



(4) 
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FIG. 4: Chemical potential versus water density for u/v = 1 and T = 0.05. The coexistence 
between a gas (p ~ 0) and a low density phase [p ~ 0.75) can be seen from the van der Waals loop. 
Note the abrupt increase for p ~ 0.75. 




Under the the physical conditions x > 0, ?/ > and r > 0, since x, y and r are ratios 
between partition functions, we must have x = y. EqsE] reduce to two equations which are 
solved numerically for specific values of the model parameters u/v. Chemical potential vs 
density isotherms are then obtained from EqEl as shown in Figures 4 and 5. 

The first striking result is that temperature destabilizes the liquid-liquid transition. Van 
der Waals loops are present only for the gas-liquid, but not for the liquid-liquid transition. 
However, for larger bond energies u/v > 0.5, the liquid phase is of low density (p = 0.75), 
whereas for smaller bond energies, the liquid phase is of high density (p = 1). A typical 
isotherm for the first case, for which the gas-liquid transition is to a phase of low density, is 
shown in Fig 4, for u/v = 1 . It is to be noted that the chemical potential presents a very 
abrupt rise, just above the gas-liquid transition, as if signaling a quasi liquid- liquid transition. 
Fig 5 shows a typical low temperature isotherm for the second case, for u/v = 0.25, for which 
no low density phase is found. 

Inspite of the absence of a liquid-liquid transition, present in the results of simulations 



see 



40| and next section), both the LD and the HD phases appear, for different strengths 



of the hydrogen bonds. Also, the abrupt rise of the chemical potential in Fig. 4 indicates 
that this mean-field treatment comes near to but misses the LD-HD transition. Similar 
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FIG. 5: Chemical potential versus water density for u/v = 0.25 for T = 0.4. The LD phase is no 
longer present, since coexistence is between a gas (p ~ 0) and a high density phase (/o ~ 1). 

problems were noted for bonding ice-type models Cyclic ordering of edges around 

a vertex, in the lattice ordered structure, cannot be represented on the tree, which may lead 
to analytical free-energies for the same parameters for which exact |45|] or simulation j4Qj| 
studies predict phase transitions. 



IV. THE MONTE CARLO RESULTS 



The model properties for finite temperatures were obtained through Monte Carlo simu- 
lations in the grand-canonical ensemble using the Metropolis algorithm. Particle insertion 
and exclusion were tested with transition probabilities given by w{insertion) = exp{~A(j)) 
and w{exclusion) = 1 if A0 > or iv{insertion) = 1 and w{exclusion) = exp{+A(f)) if 
A(f) < with A(f) = exp{i3{eparticie — 1^) — ln(18)} where Cparticie is the particle energy. Since 
the empty and full sites are visited randomly, the factor 18 is required in order to guarantee 
detailed balance. 

A detailed study of the model properties was presented in a previous publication for a 
particular value of the energy parameter, u/v = I ■4fl|. Here, we present the evolution of the 
phase diagram under variation of this energy parameter which represents relative energy of 
bond strength. Note that the relevant parameter range is u/v > 1/2. For u < v/2, the LDL 
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FIG. 6: Chemical potential vs. density isotherms for different temperatures, p is given in units of 
lattice space and the T = kBT/v. Here we illustrate the case u/v = 1. 
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FIG. 7: Phase-diagram showing chemical potential vs. temperature for u/v = 0.55,0.6,0.75,1. 
The solid lines are the LDL-HDL coexistence lines. The gas-LDL collapse into a single dashed line. 

The cocxistGnce at zero temperature at n/v = —1.6,-1.2,0,2 and at ^/v = —2 are exact. The 
error bars refer to the size of the hysteresis loop. For visualization purposes the error bars for the 
gas-liquid points are not shown. The statistical error bars have similar size as the symbols. 
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FIG. 8: Phase-diagram showing pressure vs. temperature for ujv = 0.55,0.6,0.75,1. Reduced 
pressure p/v is given in units of lattice space. The soUd hnes ( and filled symbols) are the LDL- 
HDL coexistence lines. The empty symbols indicate the gas-LDL coexistence lines collapsed into 
a single dashed line. The coexistence at zero temperature at p/v = 0.3, 0.6, 1.5, 3 and p/v = are 
exact. 

disappears. Our study is for L—10 and runs were of the order of 10^ Monte Carlo steps. 

Chemical potential vs. density isotherms for u/v — 1, and a large set of temperatures 
are illustrated in Fig.6. The figure shows the first-order phase transitions between the 
gas and the low density liquid phase and between the low density liquid and high den- 
sity liquid phases. The coexistence lines arc estimated from the mid-gaps in the chemical 
potentials. Fig. 7 illustrates the chemical potential vs. temperature phase-diagram for 
u/v — 0.55, 0.6, 0.75, 1. The gas-LDL coexistence lines for the different u/v seem to collapse 
into a single line that at zero temperature gives /i/v — —2. The different LDL-HDL lines 
are also illustrated. At zero temperature, exact calculations locate the end-points of these 
lines at, respectively, n/v — —1.6, —1.2, 0, 2. The error bars here are a measure of the size 
of the hysteresis loops for each temperature. The statistical errors are of similar size as the 
symbols. As w/f ^ 0.5 the liquid-liquid coexistence lines approach the gas-liquid transition, 
merging with it at u/v = 0.5. 

Pressure was computed from numerical integration of the Gibbs Duhem equation, at fixed 
temperature, from zero pressure at zero density. Fig. 8 shows the corresponding pressure 
vs temperature phase diagrams. 
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At zero temperature, we have the exact values p/f = 3, 1.5, 0.6, 0.3. An inversion of the 
LDL-HDL phase boundary, close to the critical point, is to be noted. The LDL-HDL line 
has positive slope for u/v = 1,0.75, and negative slope for u/v = 0.55. This indicates that 
for u/v = 1,0.75 the low density liquid phase is more entropic than the high density liquid 
phase as it would be the case for most liquids. However, for u/v = 0.55, the LDL phase 
is less entropic than the HDL phase, a behavior that is expected for water. The line of 
temperature of maximum densities, TMD, is present for all u/v > 0.5 values. The figure 
illustrates the TMD only for the u/v = 1 case. 

In summary, we have shown, both from mean-field and from simulations that the model 
we proposed in a previous publication may present two liquid phases and the associated 
double criticality, depending on the ratio of bond strength to energy penalty for distortion, 
represented by our parameter u/v. Double criticality disappears in case the price paid for 
increasing the coordination is too high. A second consequence of varying the relative energies 
is that the slope of the liquid-liquid line may become negative, around the critical point. 

Analysis of the model behavior shows that the two liquids are present if the two competing 
distances are sufficiently separated, but the "shoulder" need not be negative (see Fig. 2), 
in contradiction to a recent proposal 47[. In other words, two characteristic attractive 
distances might also generate the liquid-liquid transition. 

Experimental and numerical evidence suggests that water possesses a first-order transition 
between a low density liquid and a high density liquid. The present model shows two liquid 
phases and consequently two critical points and a line of density anomalies. 

Experimental evidence for a first-order transition in silica is still missing. According to our 
model, as the bond energy varies, the critical point temperature decreases. For temperatures 
above the critical point, under compression the liquid goes from the low density phase to 
high density phase in a continuous way as it is observed experimentally in silica. Since the 
liquid-liquid critical point moves to lower temperatures as the bond energy becomes weaker, 
one could suggest that the first-order transition predicted by simulations is not observed in 
experiments because the coexistence line and the second critical point in the case of silica 
might be located at very low temperatures. 

Water and silica are network bonding systems, and present low-coordination solids with 
transition to higher coordinated solids at high pressures. Ice I is tetrahedrally bonded, with 
bonds suffering distortion, for higher pressures. At still higher pressures, the tetrahedral 
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structures interpenetrate, to satisfy both energy and pressure requirements. Silica transitions 
between tetrahedrally and octahedrally coordinated structures. It would be interesting to 
establish energy criteria for bond distortions in both cases, in order to test the ideas we 
propose in this study. 
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